# -*- coding: utf-8 -*-
"""
Created on Wed May 30 13:53:17 2012
af optimization test
@author: Maxim
"""

import numpy as ny
import ga2
import airfoil as Airfoil_Class
import matplotlib.pyplot as plt
from datetime import datetime

def CostFcn(X):
    Au = X[0:4]
    Al = X[4:7]
    Al = ny.hstack([-Au[0],Al])
    Mach = 0.16
    Re = 4.17e6
    alphaStart = -10
    alphaEnd = 20
    alphaStep = 1.0
    clreq = 0.4

    af = Airfoil_Class.Airfoil()
    af.CST_airfoil(Au,Al)
    tc = ny.max(af.upOrig[:,1]-af.loOrig[:,1])

    g = ny.zeros([4])
    af.calc_Xpolar(Mach,Re,([alphaStart,alphaEnd,alphaStep]),50)

    g[0] = 0.14 - tc
    g[1] = tc - 0.18
    g[2] = 1.6 - af.polar.getCLmax()
    f,alpha = af.polar.cd_at_cl(clreq)
    g[3] = -alpha
    h = []
    return f, g,h

def getGradient(Fcn,x0):
    dx = 1e-3
    dimm = x0.shape[0]
    f = Fcn(x0)[0]
    grad = ny.zeros([dimm,1])
    for ii in range(dimm):
        tmp = ny.copy(x0)
        tmp[ii] += dx
        df = Fcn(tmp)[0] - f
        grad[ii] = df/dx
    return grad

StartTime = datetime.now()

bestAfpath = r'D:\codes\actools\pyAC\actools\results'
resultFile = bestAfpath + '\\result1.dat'
airfoilFile = bestAfpath + '\\airfoil1.dat'
historyFile = bestAfpath + '\\history.dat'

lb = ny.array([0.08,0.001,0.001,0.001,-0.4   ,-0.4   ,-0.4])
ub = ny.array([0.4 ,0.4  ,0.4  ,0.4  , 0.001, 0.001, 0.001])

opt = ga2.gaOptions(lb,ub)
opt.CorrNan = True
opt.PopSize = 75
opt.MaxIterations(50,30)
opt.sigmaEnd = .1
opt.sigma = 1.0
opt.HistFile = historyFile
fBest,xBest,fHist,xHist,Iter = ga2.gaMain(CostFcn, opt)

Au = xBest[0:4]
Al = xBest[4:7]
Al = ny.hstack([-Au[0],Al])
optimizedAf = Airfoil.Airfoil()
optimizedAf.CST_airfoil(Au,Al)
optimizedAf.write_airfoil_txt(r'D:\4. Conferences\ATIO - sep 2012\resutls\gaOptimum.dat')
f,g,h = CostFcn(xBest)

TimeConsumed = datetime.now()-StartTime

print fBest, xBest
print 'f:   ',f,'\ng:    ',g
print 'Time consumed:', TimeConsumed

optimizedAf.plotAirfoil()

print getGradient(CostFcn, xBest)